Full magnetoelectric response of Cr 2 03 from first principles 
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The linear magnetoelectric response of Cr2C>3 at zero temperature is calculated from first prin- 
ciples by tracking the change in magnetization under a macroscopic electric field. Both the spin 
and the orbital contributions to the induced magnetization are computed, and in each case the 
response is decomposed into lattice and electronic parts. We find that the transverse response is 
dominated by the spin-lattice and spin-electronic contributions, whose calculated values are consis- 
tent with static and optical magnetoelectric measurements. In the case of the longitudinal response, 
orbital contributions dominate over spin contributions, but the net calculated longitudinal response 
remains much smaller than the experimentally measured one at low temperatures. We also discuss 
the absolute sign of the magnetoelectric coupling in the two time-reversed magnetic domains of 
Cr 2 3 . 

PACS numbers: 75.85.+t,75.30.Cr,71.15.Rf,71.15.Mb 
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I. INTRODUCTION 

There has been a recent resurgence of interest in mag- 
netoelectric (ME) couplings in solids. 1 Of particular im- 
portance is the linear ME effect, which can occur in insu- 
lating materials with broken inversion and time-reversal 
symmetries. It can be described by a response tensor 
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where P is the electric polarization induced by the mag- 
netic field H, and conversely M is the magnetization in- 
duced by the electric field £. 

The early milestones in the long history of the linear 
ME effect include the original prediction by Dzyaloshin- 
skii that it should occur in Ci^Os, 2 and its observation 
shortly after, both in M(£) 3 ' 4 and in P(H) 5,6 measure- 
ments. The ME effect has since been observed in a large 
variety of materials, but Cr203 remains one of the best- 
studied ME compounds. The early literature is surveyed 
in the monograph by O'Dell, 7 and recent reviews are 
given in Rcfs. 1, 8-10. 

Most of the early theoretical work was phenomenologi- 
cal in character, making it difficult to assess the dominant 
mechanisms behind the ME response. These can be di- 
vided into electronic (i.e. frozen-ion) vs. lattice responses 
on the one hand, 11 and spin vs. orbital magnetic contri- 
butions on the other. 12 Ab initio theory is an ideal tool 
for unraveling the microscopic mechanisms of the ME ef- 
fect in real materials, and the first calculations started 
to appear in recent years. The initial focus was on spin- 
lattice contributions, 13 ' 14 in part because investigations 
of related phenomena in multiferroic materials over the 
last decade had indicated that spin-lattice effects are of- 
ten dominant there. 15 In reality, however, very little is 



known about the relative magnitudes of the various con- 
tributions to the ME tensor in typical magnetoelectric 
materials. 

Evidence for a significant electronic ME response in 
&2O3 came from optical measurements at frequencies 
above the lattice resonances: in a series of milestone 
experiments, 16-18 Pisarev, Krichevtsov, and collabora- 
tors observed optical effects governed by an effective ME 
tensor ot(u>), and found it to be comparable to the static 
ME coupling. Regarding the distinction between spin 
and orbital couplings (e.g., how much of the £-field in- 
duced magnetization comes from spin moments versus 
orbital currents), it is probably rather difficult to sepa- 
rate them experimentally due to the weakness of the ME 
effect in known ME materials. Investigation of the orbital 
contribution to the ME response is however interesting 
in its own right. In particular, it was recently established 
that Z2 topological insulators with broken time-reversal 
symmetry on the surface should display a quantized elec- 
tronic orbital ME response 19 ' 20 with a relatively large 
quantum (a = 24.3 ps/m in SI units). This result further 
suggests that large orbital ME responses can in principle 
be achieved even in generic (non-topological) insulators 
with strong spin-orbit coupling without any constraint 
on surface preparation. 21 

In this paper, we carry out a thorough first-principles 
investigation of the linear ME effect in the paradigmatic 
system Cr 2 03. We compute the full static response, in- 
cluding on the same footing all four basic contributions: 
spin-lattice, spin-electronic, orbital-lattice, and orbital- 
electronic. This completes the programme initiated in 
Refs. 13 and 22, where some but not all of them were 
evaluated. As in those works, we shall focus exclusively 
on the ME response at zero temperature, which is de- 
termined by mechanisms involving the spin-orbit inter- 
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FIG. 1. (Color online.) Rhombohedral primitive cell of 
Cr203. The arrows indicate the orientations of the magnetic 
moments on the Cr ions (the spins on the ions point opposite 
to the arrows). The center of the cell is a center of inversion 
symmetry coupled with time-reversal. 

action. 

We find that for the response transverse to the rhom- 
bohedral axis the spin contributions are much larger than 
the orbital ones. The calculated values of the lattice and 
electronic responses are in good agreement with both 
static and optical ME measurements, as well as with 
previous calculations. In the case of the longitudinal 
response the calculated orbital contributions are larger 
than their spin counterparts in both the electronic and 
lattice channels. However, as a result of a near cancel- 
lation between the orbital-electronic and orbital-lattice 
contributions, the total calculated longitudinal response 
is negligibly small. Thus, the nonzero longitudinal re- 
sponse that is measured at low temperatures remains 
unaccounted for. Some possible reasons for this disagree- 
ment will be discussed. 



II. PRELIMINARIES 

A. Structure of Cr203 

Chromium (III) oxide (eskolaite) crystallizes in a 
corundum- type structure shown in Fig. 1, with two for- 
mula units per primitive cell. The magnetic space group 
is R3cl' above TV = 307 K. Below this temperature 
Ci'203 turns into an antiferromagnetic (AFM) insula- 
tor, with magnetic space group R3V. The magnetic 
moments on the Cr ions align along the rhombohedral 
z axis, pointing up and down in an alternating manner 
(see Fig. 1). The magnetic point group is 3 to', which 
allows for a diagonal ME tensor a. with two independent 
components, a±_ = a xx = a yy (transverse) and ay = a zz 
(longitudinal). 23 

We note that there are two distinct possibilities for ar- 
ranging the magnetic moments in the AFM ground state, 
related to one another either by time reversal (i.e., by flip- 
ping the magnetic moments on every Cr ion) or by spa- 



tial inversion. As each of these operations also flips the 
sign of a, it is important to specify which configuration 
is assumed when reporting values for a± and an. Our 
calculations refer to the configuration shown in Fig. 1. 



B. Formalism and review of previous calculations 

We begin by clarifying issues of units and conventions. 
Equation (1), which is written in the (£,H) frame, con- 
forms with the standard experimental definition of the 
linear ME tensor, which has units of ps/m in SI units. 
Instead, from the point of view of first-principles theory 
it is more convenient to work in the (£ , B) frame, where 
a has units of vacuum admittance \Ao//zo, 

( dPi\ fdMj\ 

The two definitions Eqs. (1) and (2) are related by 
ot £H = fj,a, where fi is the magnetic permeability. In the 
approximation that fi/fio — 1, which is a good approxi- 
mation for most non-ferromagnetic materials, the conver- 
sion is trivial, and we shall report the calculated values of 
a as though we had computed them in the (£, H) frame. 
For a more detailed discussion, see Sec. II. A of Ref. 21. 

Let us now discuss how to compute the various con- 
tributions to the ME tensor. To recap, the full response 
can be decomposed into spin and orbital parts according 
to the nature of the induced magnetization in the M(£) 
picture. Each of these can be further decomposed ac- 
cording to the two basic mechanisms by which the field 
acts on the system. The electronic part describes the 
ME response that the system would have if the ions were 
held fixed in their equilibrium positions. The remain- 
ing lattice part is associated with the field-induced ionic 
displacements. 



1. Lattice response 

We consider first the calculation of lattice couplings. 
Here the influence of the applied field (B or £) on the 
non-conjugate moment (P or M) is mediated by internal 
ionic displacements u, so that 

latt = apau = /9M9u\ T 

dudB \ du d£J ' 1 ' 

where a summation over the atoms in one crystal cell is 
implied, and 'T' denotes the matrix transpose. (In gen- 
eral there may also be a strain mediated coupling, 14 but 
in Cr203 this contribution vanishes by symmetry, and 
it will not be considered further here.) Optionally, one 
may also take advantage of the fact that the displace- 
ments induced by the field are mediated by field-induced 
forces F. For the case of applied electric field, Eq. (3) 
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can be rewritten as 13 
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where E is the total energy per unit cell and fi is the 
unit-cell volume. Here we have made use of the fact that 
the Born effective charge tensor can be expressed equiva- 
lently as (dF/d£) T = fldP/du. [Alternatively, by invok- 
ing the magnetic analog (9F/<9B) T = fldM./du we can 
arrive at this same equation in a different way, starting 
from Eq. (3) for the case of applied magnetic induction.] 
Note that the inverse of the force-constant matrix now 
appears symmetrically between the magnetic and electric 
Born tensors in Eq. (4). 

There are several choices on how to proceed. One pos- 
sibility, following Eq. (3), is to relax the structure in 
the presence of a small B or £ field, and then compute 
the relaxation-induced change in P or M. Alternatively, 
Eq. (4) expresses ct latt in terms of three basic quanti- 
ties (the force-constant matrix, the Born charges, and 
their magnetic analogs), all of which can be computed 
as changes of various quantities in response to atomic 
displacements at vanishing fields. One can choose to 
compute such derivatives by finite differences or by us- 
ing linear-response techniques available in most density- 
functional packages. 



2. Electronic response 

The calculation of the electronic response ct cl requires 
coupling the field B or £ in Eq. (2) directly to the elec- 
trons, and determining the induced P or M. In practice 
this can be done using either finite-field approaches or 
linear-response techniques. 

Of the two contributions, spin-electronic and orbital- 
electronic, the latter is the most challenging one to calcu- 
late. A perturbative expression valid for periodic crystals 
was recently derived, 24 ' 25 which can be implemented in 
the context of density-functional perturbation theory. In 
the present work we have opted to calculate the orbital- 
electronic response as <9M orb /9£, using finite electric 
fields. Another possibility would be to calculate it as 
<9P/<9B orb , using a finite orbital magnetic field. The in- 
clusion of orbital magnetic fields in total-energy calcula- 
tions of periodic solids is, however, a challenging problem 
which has not yet been fully solved, in spite of some re- 
cent progress 
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3. Review of previous calculations for CV-2 O3 

The methods described above were recently used to 
evaluate the spin-lattice and spin-electronic parts of a. 
For the spin-lattice contribution, Ihiguez 13 performed his 
pioneering calculations following Eq. (4), while Bousquct 
et al. 22 used Eq. (3). More precisely, the latter authors 



performed structural relaxations in the presence of a fixed 
Zeeman magnetic field B spm by adding to the Kohn- 
Sham energy functional a Zeeman term describing the 
coupling to the spins. Furthermore, by monitoring the 
linear change in the electronic polarization P cl under a 
small field with fixed ions, Bousquet et al. were also 
able to determine the spin-electronic response. Thus, 
out of the four possible contributions to a, only the two 
spin contributions (lattice and electronic) have previously 
been evaluated from first principles for C^Oa. 



C. Computational approach 

Let us now describe the method that we use for calcu- 
lating the lattice and electronic ME responses, including 
in each case both the spin and the orbital parts of the 
response. 

For the lattice couplings we employ a method similar 
to that of Ref. 13 but including also the orbital contri- 
bution to <9M/<9u (we found this to be a more efficient 
approach than relaxing the lattice under a finite electric 
field). We first compute the Born effective charges and 
force-constant matrix using linear-response techniques, 28 
and from these we find the first-order field-induced dis- 
placements Au = (du/d£) ■ A£, where a nominal field 
A£ of ~ 10 9 V/m is applied along the rhombohcdral axis 
or in the perpendicular direction. Displacing the atoms 
by Au, we then determine the induced magnetization 
AM s P in + AM orb . The linearity of the magnetization re- 
sponse was checked by both reducing the magnitude and 
flipping the sign of A£. 

In order to reduce the computational cost, the spin- 
orbit interaction is not included in the linear-response 
calculations. This procedure captures the dominant con- 
tributions to a latt , i.e., those that are linear in the spin- 
orbit coupling strength; we have checked that it pro- 
duces results which are almost identical to a calculation 
in which the spin-orbit coupling is included at every step. 

For a given set of ionic displacements Au, the or- 
bital magnetization at £ = is calculated under periodic 
boundary conditions as 



where 29 31 
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LC 



M orb = M LC + M IC , 
d 3 k 



(5) 
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Here 'LC and 'IO stand for local circulation and itin- 
erant circulation, respectively, |w„k) is the cell-periodic 
part of the Bloch state IV'nk), and = e~ ik r %e lk ' r , 
where % is the Kohn-Sham Hamiltonian of the crys- 
tal. Summations over occupied bands are implied for 
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repeated band indices, and Vk = (1 — |u„k)(itnk|)Vk. 
In practice the Brillouin-zone integral is replaced by a 
summation over a uniform grid, and Vk is evaluated 
on that grid by finite differences. 31 M LC and M IC are 
separately gauge-invariant, i.e., they remain unchanged 
under fc-dependent unitary transformations among the 
occupied states. 

Let us now turn to the electronic response, which 
we calculate as 9M/<9£, taking advantage of the well- 
established ab initio treatment of homogeneous elec- 
tric fields in periodic insulators. 32 The magnetization 
M spm + M orb is determined with and without an elec- 
tric field of intensity ~ 10 9 V/m (using in both cases 
the same crystal structure optimized at zero field) in or- 
der to extract the spin-electronic and orbital-electronic 
ME couplings. The evaluation of spin magnetization is 
straightforward and here we just mention that, as an ad- 
ditional check, we have recomputed the spin-electronic 
coupling using the converse Zeeman-ficld approach, find- 
ing good agreement between the two methods. 

To compute M orb at finite £, we make use of the fol- 
lowing generalization of Eqs. (5-7). 24 One part is given 
by the same expression valid at zero field, Eqs. (6) and 
(7), upon reinterpreting the states \u n k) therein as field- 
polarized Bloch states 32 (and ~K as the crystal Hamilto- 
nian calculated from the field-polarized periodic charge 
density). To this, an additional contribution of the form 
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must be added in order to obtain the full orbital magne- 
tization. Here A£ m = i(u n ^\ Vk|u n k) is the Berry con- 
nection matrix; the integrand is a scalar known as the 
Chcrn-Simons 3- form 19,20 (band indices are suppressed). 
Thus, at £ ^ we have, instead of Eq. (5), 



M 



orb 



M 



LC 



M 



ic 



M 



cs 



(9) 



and all three terms contribute to the orbital-electronic 
ME response. The term a cs is purely isotropic and 
can be calculated from the valence Bloch states at zero 
field. Its numerical evaluation requires a smooth gauge in 
k space, and this can be achieved by mapping the valence 
bands onto localized Wannicr functions. 21 



D. Technical details 

The total-energy and linear-response calculations were 
performed using the Quantum-ESPRESSO 33 ab initio 
code package, working in a fully relativistic framework 
where the spin-orbit interaction is included in the atomic 
pseudopotcntials. We employed Troullicr-Martins norm- 
conserving pseudopotentials, 34 which in the case of Cr 
included the semi-core 3s and 3p states in the valence. 

The wavefunctions in the solid were expanded in plane 
waves with an energy cutoff of 250 Ry for structural relax- 



TABLE I. Calculated and experimental structural parame- 
ters of Cr2C>3 in the antiferromagnetic phase: rhombohedral 
lattice parameter a, rhombohedral angle a, and Wyckoff po- 
sitions of the Cr ions (4c orbit) and O ions (6e orbit). 





a (A) 


a (deg) 


Wyckoff positions 








Cr 


O 


PBE (This work) 


5.415 


54.45 


0.1541 


0.0597 


LDA (Ref. 21) 


5.322 


53.01 


0.1575 


0.0690 


Expt. (Ref. 35) 


5.358 


55.0 


0.1528 


0.0566 



ations and linear-response calculations and 150 Ry for or- 
bital magnetization calculations. The Brillouin zone was 
sampled on a 4 x 4 x 4 Monkhorst-Pack mesh for most self- 
consistent-field (SCF) calculations. While this mesh den- 
sity produced converged values for the spin-lattice and 
spin-electronic ME contributions, the two orbital con- 
tributions converged more slowly with fc-point sampling 
[this is probably related to the finite-differences represen- 
tation of the covariant derivatives in Eqs. (6) and (7)]. 
After testing several grid densities, we concluded that a 
7x7x7 mesh gave sufficiently converged values. 



As noted in Ref. 13, the computation of ME couplings 
demands a very tight tolerance on the convergence of 
the self-consistent field loop. We therefore used rather 
stringent convergence thresholds, of the order of 10~ n - 
10~ 12 Ry in the total energy. In order to reach this 
level of convergence in a reasonable number of steps with 
Quantum-ESPRESSO, we found it useful to use the 
Thomas-Fermi charge mixing scheme, 36 by setting the 
input variable 'mixingjnode ' to 'local-TF'. We also 
found that the speed of convergence of the calculations 
with a finite electric field was improved by increasing the 
field gradually from zero in small steps. 



The exchange-correlation potential was described 
within the generalized-gradient approximation 
(GGA) using the Perdew-Burke-Ernzerhof (PBE) 
parametrization. 37 This choice was made after having 
optimized the structure using both the local-density 
approximation (LDA) and PBE, and finding that the 
latter produced structural parameters in better agree- 
ment with experiment (see Table I). In particular, LDA 
underestimates the unit-cell volume by 7.3% while PBE 
overestimates it by only 1.7%. We note that the authors 
of Refs. 13 and 22 used LDA+U with the experimental 
cell volume enforced. As for the magnetic structure, 
the staggered spin moments on the Cr atoms have a 
value of 2.7/ZB/atom, for a sphere integration radius of 
1.3 A. This is in good agreement with the LDA+U value 
reported in Ref. 13. 
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TABLE II. Calculated contributions to the magnetoelectric 
tensor components a± and Qy in Cr 2 03 . Columns (rows) 
show the spin and orbital (electronic and lattice) contribu- 
tions. [The results from previous calculations are indicated in 
parentheses.] 



c«i (ps/m) 
Spin Orb. Total 



a || (ps/m) 
Spin Orb. Total 



Elec. 



0.26 -0.014 0.25 



Latt. 



(0.34 a ) 

0.77 
(l.ll a ) 
(0.43 b ) 
Total 1.03 



0.025 0.80 



0.011 1.04 



0.0007 -0.009 -0.008 

(o a ) 

0.0026 0.008 0.011 

(o a ) 

(0.00 b ) 

0.003 -0.001 0.002 



a Ref. 22. 
b Ref. 13. 



III. RESULTS 
A. Contributions to the ME response 

The main results of our calculations are presented in 
Table II together with results from previous theoretical 
works, given in parentheses. Let us first analyze the 
transverse ME response. The magnitude of the calcu- 
lated static value, |a_i_| = 1.04 ps/m, agrees well with 
the most reliable measurements, which range from 0.7 to 
1.6 ps/m. 38 ' 39 The spin-lattice contribution accounts for 
about 75% of that value, with the remaining 25% coming 
mostly from the spin-electronic response, while the two 
orbital contributions are negligible (less than 2%). The 
values we obtain for the individual contributions a 1 ^ 1 
and af_ agree well with those calculated in Ref. 22 using 
the converse Zeeman-field approach. 

In the case of the longitudinal response, the relative 
strengths of the four contributions are very different. As 
in previous calculations, 13 ' 22 we find that the spin contri- 
butions to a || are very small, summing to only 0.003 ps/m 
in our calculation. This can be understood as resulting 
from the extreme stiffness of the magnitude of the spin 
moment in a collinear band antifcrromagnet, which is 
also reflected in the near- vanishing of the spin magnetic 
susceptibility %|| at T = 0. 12 

Experimentally, however, the low-temperature a\\ is 
found to be about 0.2-0.3 ps/m. 38 ' 39 This is smaller than 
a± by a factor of 3 to 6, but still about two orders of mag- 
nitude larger than our theoretical spin value, suggesting 
that orbital effects might be responsible for most of the 
a|| response. Indeed, Hornreich and Shtrikman 12 pointed 
out that a zero-temperature longitudinal ME response 
could arise in Cr2C>3 from an electric-field- induced shift 
in the g factor of the Cr ions (see also Ref. 11). This is 
an orbital effect that should be automatically included in 
the present calculations. In fact, we do find that our com- 
puted orbital-lattice and orbital-electronic contributions 



to au are nearly an order of magnitude larger than the 
corresponding spin contributions. However, the orbital- 
lattice and orbital-electronic contributions individually 
are still an order of magnitude smaller than the mea- 
sured value. Moreover, these two contributions have op- 
posite signs, resulting in a near cancellation of the entire 
longitudinal response. Our total ay of 0.002 ps/m thus 
remains about two orders of magnitude smaller than the 
measured value. 

There are several possible explanations for this discrep- 
ancy. First, the theoretical values for the orbital longi- 
tudinal response are quite small, and thus might be es- 
pecially sensitive to numerical errors. However, we have 
checked /c-point and self-consistent convergence carefully, 
and do not believe this is a major concern. More seri- 
ous is the potential dependence on choice of exchange- 
correlation potential. In particular, within the LDA we 
found that the orbital-electronic and orbital-lattice con- 
tributions are approximately a factor of 3 larger com- 
pared to PBE, although similar cancellation of the two 
contributions was observed. Future work is needed to 
check the sensitivity of these calculations to the choice 
of GGA (adopted here) as opposed to LDA, LDA+U or 
GGA+U, hybrid functional, or other orbital-dependent 
functionals. Since orbital currents play a crucial role, the 
use of current-density functionals should probably also be 
explored. On the experimental side, it would probably be 
advisable to check the dependence of the measured value 
on sample quality, in order to rule out extrinsic effects 
associated with defects, surfaces, contacts, etc. 

It is also possible, however, that the experimentally 
observed response is dominated by some physics not cap- 
tured by LDA or GGA approximations to the exact den- 
sity functional. For example, the strong dependence of 
a || upon temperature makes it clear that thermal fluctu- 
ations strongly influence the longitudinal response. By 
the same token, it is possible that quantum spin fluctu- 
ations, already present in the antiferromagnetic state at 
zero temperature, may play an important role. For the 
time being, we leave this as an open question. 

Before closing this section, we recall that the orbital 
ME response can be further decomposed into local circu- 
lation (LC), itinerant circulation (IC), and - in the case of 
the orbital-electronic response - Chern-Simons (CS) con- 
tributions, as in Eq. (9). Table III shows the breakdown 
of the full orbital response computed in the present work. 
In our previous study of Cr203, only the isotropic CS 
term was calculated (using LDA rather than GGA). 21 In 
that work we found the CS term to be <~ 0.01 ps/m, about 
an order of magnitude larger than the presently calcu- 
lated value. Further work is needed to determine how 
the various terms in the ME response of Cr 2 03 depend 
on the choice of exchange-correlation potential. It can be 
seen that the CS contribution to the orbital-electronic re- 
sponse is about an order of magnitude smaller than the 
LC and IC contributions. Individually, the LC and IC 
orbital-electronic contributions are somewhat larger for 
ev|| than for a±, but taken together the opposite is true. 
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TABLE III. Decomposition of the calculated orbital ME re- 
sponse of Cr 2 03 (presented in Table II) into "local circula- 
tion", "itinerant circulation", and "Chern-Simons" contribu- 
tions coming respectively from Eqs. (6), (7), and (8). 



a° Ib (ps/m) 



Electronic 

Local circulation 

Itinerant circulation 

Chern-Simons 

Subtotal 
Lattice 

Local circulation 

Itinerant circulation 

Subtotal 
Total 



-0.0064 
-0.0084 
0.0012 
-0.0136 

0.0202 
0.0051 
0.0253 
0.0117 



-0.0237 
0.0135 
0.0012 

-0.0090 

0.0078 
0.0000 
0.0078 
-0.0012 



As for the orbital-lattice contributions to a± and au , they 
come mainly from the LC terms. 



B. Sign of the ME response 

We now discuss the overall sign of the tensor a. As 
already mentioned, in Cr 2 03 this sign depends on the 
orientation of the magnetic moments (see Fig. 1). Ex- 
perimentally, a single AFM domain can be stabilized by 
cooling the sample through the Neel temperature in the 
presence of parallel (or antiparallel) electric and mag- 
netic fields ( "magnetoelectric annealing"), and the spin 
structure can then be analyzed using spherical neutron 
polarimetry. 40 ' 41 

According to Ref. 40, the orientation of the magnetic 
moments shown in Fig. 1 therein corresponds to a do- 
main annealed with electric and magnetic fields pointing 
in the opposite direction along the rhombohedral axis, 
provided that arrows in that figure indeed indicate direc- 
tions of spin moments rather than magnetizations. Since 
the magnetoelectric tensor appears in the free energy in 
the form Fme = —(XijSiHj, the domain under considera- 
tion should have negative ot\\ near the Neel temperature. 
Experimental measurements of magnetoelectric coupling 



as a function of temperature 38,39 show that a\\ changes 
sign around 100 K, while a± is negative all the way to 
4.2 K. Assuming that magnetic domain is determined at 
high temperatures, close to the Neel temperature, and 
that magnetic domains remain frozen upon cooling to 
4.2 K, we can conclude that at 4.2 K the domain shown 
in Fig. 1 must have a± > and au > 0. 

Our computed signs appear to agree with the experi- 
mental work of Ref. 40, although it was not made entirely 
clear whether the signs reported there refer to spins or 
magnetizations. Now that first-principles theory is seri- 
ously beginning to confront experiment in the field of 
magnetoelectric couplings, we urge closer attention to 



sign issues in future investigations, both theoretical and 
experimental. 



C. Comparison to optical measurements 

We now turn to the comparison with existing mea- 
surements of the optical ME tensor a(uj). As our theory 
only deals with static fields, the calculated a el should be 
thought of as the uj — > limit of the purely electronic 
optical response (quasistatic limit). This is expected to 
approximate reasonably well the measured response at 
frequencies between the lattice and electronic resonances 
and sufficiently far from both. 

The ME coupling influences both the transmission and 
reflection of light from a magnetoelectric medium, giving 
rise to characteristic optical effects which are odd under 
time reversal. 11 ' 42 While the propagation of electromag- 
netic waves inside a ME medium is only affected by the 
traceless part of a, all tensor components can in principle 
be extracted from reflectance measurements, although in 
that case the net effect may also have surface-specific 
contributions. 17 The reflection experiments of Ref. 17 
were carried out using visible light with a wavelength 
of 633 nm (1.96 eV), which falls within the exciton ab- 
sorption range of Cr203, thus precluding a meaningful 
comparison with our quasistatic calculations. 

We therefore focus on the earlier transmission 
measurements, 16 which used near infrared light of 
1156 nm (1.07 eV). The effect that was observed consists 
of a tilt away from the crystallographic y and z direc- 
tions of the linear polarization of light traveling along x. 
The tilt angle <f) is related to the components of the op- 
tical ME tensor (expressed in Gaussian units) and index 
of refraction by 17 ' 42 



1 a. 



2 n z 



(10) 



While an effect which changed sign between time- 
reversed samples was clearly observed, a time-even back- 
ground signal of comparable magnitude could not be 
eliminated. The most reliable value, <j> = 4' ~ 1.2 x 
10~ 3 rad, was measured at 220-240 K. As the absolute 
value of the linear birefringence was not reported, we use 
the value n z — n x — 5.8 x 10~ 2 quoted in Ref. 17 for 
633 nm, to arrive at a xx — a zz ~ 0.12 ps/m. The agree- 
ment with our calculated value of 0.26 ps/m is quite sat- 
isfactory, given the experimental uncertainties as well as 
the limitations in our theory (namely, the DFT underes- 
timation of the optical gap and the assumed quasistatic 
and low temperature limits in the calculation). 

We emphasize that, as in the case of the static measure- 
ments discussed earlier, the dominant type of AFM do- 
main present in the samples was not specified in Ref. 16. 
Hence the absolute sign of the measured optical ME co- 
efficient was not determined. It would be interesting to 
carry out optical and static ME measurements on the 
same single-domain sample at low temperatures. This 



7 



would allow one to extract the relative sign between aj_ 
and which we predict to be positive. 

IV. CONCLUSIONS 

In summary, we have performed a thorough investiga- 
tion of the zero-temperature ME response in Cr 2 03 us- 
ing first-principles calculations. We analyzed the lattice 
and electronic parts of the response including both spin 
and orbital magnetization contributions, being careful to 
treat all four contributions on an equal footing. In par- 
ticular, we treated the orbital response using the mod- 
ern Berry-phase theory, without introducing muffin-tin 
approximations, in which orbital currents are computed 
inside spheres around atoms. 

We have then compared the calculated values with 
static and optical measurements. Previous calculations, 
which focused on the spin contributions, had found an 



essentially null value for qn, in disagreement with ex- 
periment. We therefore set out to check whether orbital 
effects could account for the observed low-temperature 
longitudinal response, as had been proposed early on in 
the literature. Our results suggest that this is not the 
the calculated orbital responses are very small, 
consistent with a scenario of strongly quenched orbital 
moments. We hope that the present findings will stimu- 
late further investigations, both on the experimental and 
theoretical sides. 

Recently we became aware of concurrent first- 
principles studies of the orbital ME response in Cr203 43 
and LiFeP04 44 using the approximation of integrating 
orbital currents within atom-centered spheres. 

This work was supported by NSF Grant DMR-10- 
05838 and by Grant MAT2012-33720 from the Spanish 
Ministerio de Economia y Competitividad. We would 
like to thank Michael Fechner, Manish Jain, and Georgy 
Samsonidzc for useful discussions. 
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